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Abstract 

Recently, research that focuses on the rigorous understanding of the relation be- 
tween simulation and/or exact models on graphs and approximate counterparts has 
gained lots of momentum. This includes revisiting the performance of classic pairwise 
models with closures at the level of pairs and/or triples as well as effective-degree- type 
models and those based on the probability generating function formalism. In this pa- 
per, for a fully connected graph and the simple SIS (susceptible-infected-susceptible) 
epidemic model, a novel closure is introduced. This is done via using the equations 
for the moments of the distribution describing the number of infecteds at all times 
combined with the empirical observations that this is well described/approximated by 
a binomial distribution with time dependent parameters. This assumption allows us 
to express higher order moments in terms of lower order ones and this leads to a new 
closure. The significant feature of the new closure is that the difference of the exact 
system, given by the Kolmogorov equations, from the solution of the newly defined 
approximate system is of order This is in contrast with the 0{1/N) difference 

corresponding to the approximate system obtained via the classic triple closure. 
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1 Introduction: exact stochastic models on and 
of graphs 

In this Section we present two important examples that motivate our investigations and a 
possible extension where progress can be made following some of the methods introduced 
in this paper. Here, we use a dynamical system type approach, where the Kolmogorov 
equations are simply considered as a system of linear ODEs with a transition rate matrix 
with specific properties such as special tri-diagonal structure and/or well defined functional 
form for the transmission rates. For example, consider a Markov chain with finite state space 
{0, 1, . . . , N} and denote by Pk{t) the probability that the system is in state k at time t (with 
a given initial state that is not specified at the moment). Assuming that starting from state 
k the system can move to either state /c — 1 or to state k + 1, the Kolmogorov equations of 
the Markov chain take the form 

Pk = ak^ipk-i - {a-k + Ck)pk + Ck+iPk+i, k = 0,...,N. (KE) 

The first motivation for our study comes from epidemiology where a paradigm disease 
transmission model is the simple susceptible-infected -susceptible (SIS) model on a com- 
pletely connected graph with nodes, i.e. all individuals are connected to each other. From 
the disease dynamic viewpoint, each individual is either susceptible (S) or infected (/) - the 
susceptible ones can be infected at rate r if connected to an infected node, and the infected 
ones can recover at rate 7 and become susceptible again. It is known that in this case the 
2 ^-dimensional system of Kolmogorov equations can be lumped to a + 1-dimensional 
system, see [12]. 

The lumped Kolmogorov equations take again the form (IKEp with 



ak = Tk{N — k), Ck = 'jk for k = 0,...,N with a_i = ctv+i = 0. (1) 

For homogeneous random graphs, where every node has n{< N) links to other nodes in 
the network, a similar system can be written down. In this case the transition rates are as 
follows, 

N - k 

Qk = '^'^kj^ — ^5 Ck = 'yk for k = 0, . . . , N with a_i = cn+i = 0. (2) 

This equation is not exact in the sense that it is only an approximation to the exact process 
unfolding on a homogeneous random graph. This is not of immediate relevance as in this 
paper we concentrate on approximations of the Kolmogorov equations via low-dimensional 
ODEs, in particular we aim to derive approximate equations for the moments of the stochastic 
process described by the Kolmogorov equations. 

A final example stems from a recent model of a dynamic network where links can be 
activated and deleted at random while subjected to a global constraint on the total number 
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of links in the network (see Kiss et al. [7]). This example is a prerequisite to models where 
dynamics on the network is coupled with the dynamic of the network Gross and Blasius [3] . 
In the case of epidemic propagation on networks it is straightforward to assume that the 
propagation of the epidemic has an effect on the structure of the network. For example, 
susceptible individuals try to cut their links in order to minimise their exposure to infection. 
This leads to a change in network structure which in turn impacts on how the epidemic 
spreads. The first step in modelling this phenomenon is an appropriate dynamic network 
model such as the recently proposed globally-constrained Random Link Activation-Deletion 
(RLAD) model. This can be described in terms of Kolmogorov equations as follows, 

Pkit) = a{N - {k - 1)) ( 1 - ^] p.^^it) 

Pk{t) +uj{k + l)pk+i{t), 



a(N-k) ( 1-— ^ ) +ujk 

\ ' \ J^max J 



k = N. 



where Pk{t) denotes the probability that at time t there are k active links in the network 
with denoting the total number of potential edges. It is assumed that non-active links are 
activated independently at random at rate a and that existing links are broken independently 
at random at rate u. Furthermore, the link creation is globally constrained by introducing 
a carrying capacity K^"-^, that is the network can only support a certain number of edges 
as given by K'^"'^. 

Using the above notation, here 

ttk = a{N -k) ( 1 - \ , Ck = ujk, for k = 0,...,N with a_i = c^+i = 0. (3) 



The common ingredient of all the models above is the set of Kolmogorov equations. Solving 
these even numerically, more often than not, is challenging and impossible simply due to 
the large number of equations. This has led to various approaches that set out to derive 
approximate models that are low-dimensional and can capture the exact dynamics in terms 
of the expected values of some well-defined quantities. These approaches could be broadly 
classified into heuristic and rigorous or semi-rigorous. For example, the original derivation 
of pairwise models was based on heuristic arguments [B]. More rigorous approaches include 
the lumping technique that exploits system symmetries and allows a significant reduction of 
the state space \12\. Additional models such as the effective degree models [Jj |9J and the 
probabihty generating function approach by Volz [101 HI] led to excellent agreement with 
simulations without the explicit proof of convergence but with strong probabilistic arguments 
based on results by Kurtz [8J. 

In this paper, starting from the Kolmogorov equations, given by the simple SIS epidemic 
model on a fully connected graph, the evolution equations for the moments are derived and 
are interpreted in terms of and compared to the classic pairwise equations. The equations 
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for the moments are not self-contained, and a new/novel closure is proposed. This is based 
on the fact that the number of infecteds at all times is well approximated by a binomial 
distribution with time dependent parameters. The performance of the new closure is inves- 
tigated numerically in the closing Section and we show that this is superior to the classic 
closure at the level of triples used for pairwise models. 



2 Pairwise models: closure relations and their perfor- 
mance 

It is of great interest to understand how and when the full set of Kolmogorov equations 
can be approximated via low- dimensional ODEs and also to assess the performance of the 
approximate models by, for example, working out the rate of convergence of the exact/full 
model towards the solution of the approximate or mean-field model or simply estimating 
the absolute difference between the two. Pairwise models have been widely used as an 
alternative or as a companion to simulations of epidemic models on mainly homogeneous 
random graphs. While originally, the pairwise equations have been heuristically defined, 
more recently, Simon et al. [12] and Taylor et al. [13] have shown that these can be derived 
directly from the exact Kolmogorov equations and that these pairwise equations are exact 
before closure on arbitrary graphs. Focusing on the simple SIS type model, the first moment 
of the distribution is given by 

[i] = r[SI]-^[Il (4) 

where [I](t) = Yl!k=o^Pk{t)- This is not a closed equation since [SI] itself is a variable and 
an equation for this is needed. However, we can look for an approximation whereby the 
expected number of edges [SI] ~ e[s/]([/]) = [/](A^ — [/]), that is, the expected number of 
{SI) pairs is estimated in terms of the number of the expected number of infecteds [/]. This 
now leads to a self-contained equation in terms of a new approximate variable I given by 

J = Ti{N - /) - 7/. (5) 

This is the well known compartmental model for the SIS epidemic. 

The same argument can repeated by using a closure at the level of triples rather then 
pairs. In this case, the exact pairwise equations are given by 

d[I] 



dt 
d[SI] 

dt 
d[II] 

dt 

d[SS] 

dt 



r[SI]-^[Il (6) 

7([//] - [SI]) + t{[SSI] - [ISI] - [SI]), (7) 

-27[//] + 2r([/S/] + [5/]), (8) 

27 [^/] - 2t[SSI]. (9) 
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Using the well know closure given by [Ai?C] = ^zf ^^^^^p^ (see [HE]) leads to the following 
approximate system 

7 = r57-7J, (10) 



SI = -f{II~SI) + T{SSI-ISI-SI), (11) 



// = -2^11 + 2r{ISI + SI), (12) 



SS = 2^SI-2tSSI. (13) 

The focus now shifts from the derivation of the approximate model to whether and how well 
these agree with output from the exact system. More precisely, we will simply consider the 
difference between the exact solution [/] and the approximate solutions / and I in terms of 
the magnitude of \[I]{t) — I{t)\ and |[/](t) — I(t)\, and how these depend on population or 
network size. Numerical investigation reveals that the difference is of order 1/A^ for closures 
both at the level of pairs and triple. This is shown in Fig. 1, where both the exact (Eq. f lKEp 
with coefficients given by Eq. ([T])) and the approximate systems (i.e. with pairwise closure 
given by Eq. ([5]) and triple closure given by Eqs. (|T0|) -( 1T3|) ) have been solved numerically. 
This is somewhat surprising given that, at least intuitively, it is expected that the closure 
at the level of triples will be superior to that at the level of pairs. This numerical result can 
be made more rigorous at least for the closure at the level of pairs as shown by [2], |TT] . 



3 Interpreting pairwise equations and closures in terms 
of moments 

3.1 Equations for the moments 

In this Section the analysis will focus on the derivation of the equations for the moments 
and the interpretation of the pairwise equations in terms of the moments. Let us define the 
jth moment associated with the stochastic process as follows 

fe=0 ^ ^ k=0 

where N^yj = Yj with j = 1, 2, . . .. It is straightforward to derive evolution equations for the 
moments. For example, the derivative of the first moment, and in a similar way all others, 
can be given in function of other moments upon using the Kolmogorov equations flKEp . The 
derivation for the first moment is outlined below, 



N N 

= ^kPk = ^ k{ak-iPk-i - (ofc + Ck)pk + Ck+iPk+i) 

k=0 k=0 
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N 



^{kak-iPk^i - kakPk + kckPk + kck+iPk+i] 



k=0 



By changing the indices of the summation, plugging in the corresponding expressions for 
the transition rates ak and Ck (Eq. ([1])) and taking into account that a_i = cn+i = the 
following expression holds, 



TV 



Yi{t) = + ^')(^ ^^'(^ + - ^h)Pi- (15) 

i=0 

Based on our notations (see, Eq. f|T^ ). the equation above reduces to 

Fi(t) = tNYi - TY2 - 7^1 = (/3 - l)Yi - ^F2, (16) 

where /3 = tN is the linking relation between mean-field-type and network models. Using 
a similar procedure, the equation for the second moment Y2 can be easily computed and is 
given by 

Y, = 2{(3 - 7)^2 - 2^Ys + {(3 + ^)Y, - ^Y^. (17) 
Equations ( fT6|) & f|T7|) can be recast in terms of the density dependent moments yjs to give 

yi = (/3 - 7)2/1 - ?/2, (18) 

2/2 = 2(/3 - 7)^2- 2^2/3 + ^((/3 + 7)2/1-/32/2). (19) 

These equations will play a role that is similar to that of the pairwise equations, and similarly, 
these are also exact before a closure is applied, at some arbitrarily chosen moment. The above 
equations are not closed or self-contained since the second moment depends on the third and 
an equation for this is also needed. It is easy to see that this dependence of the moments 
on higher moments leads to an infinite but countable number of equations, see p]. Hence, a 
closure is needed and in the next Section we will identify how the classic closure translates 
to a closure in moments, namely expressing 1^3 as a function of Yi and Y2. 



3.2 The equivalence between pairwise and moment equations 

To be able to link the pairwise approach to the moment approach it is necessary to count 
the expected value of the singles, pairs and triples in terms of the moments. For a fully 
connected graph this is straightforward. For example, based on Eq. f|T^ . the expected 
number of infecteds is given by 



N 

[I]{t) = Y,kpkit)=Ny^{t). 

k=0 



(20) 
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Similarly, it is easy to show that similar identities for pairs and triples can be derived. For 
example, 

N ^ k k 

[Sim = Y.^{N- k)p,{t) =N'J2j;^i^- = N'{yi{t) - y2{t)), (21) 

where, A;(A^ — k) simply denotes the number of (SI) pairs on a fully connected graph with 
N nodes and k infected individuals. For triples, the calculations are equally intuitive. For 
example, the expected number of [SSI] triples can be counted by averaging over {N — k){N — 
k — l)k - the number of [SSI] triples in the presence of k infected nodes. Hence, the following 
relation holds, 

N 



[SSI]{t) = Y,{N-k){N-k-l)kpk{t) 

I \ k 

-rPk{t) 



k=0 

N 



k=0 ^ ^ ^ 

= iV^ (^(^1 - 1^ y,it) + ( - - 2 ) ym+ysit) ] . (22) 
Following the same simple procedure as above, the following relations hold 



k 




N " 






- 2 



[S]{t) = N{l-y,{t)), (23) 
[Iim = (^y,{t) - j^y,{t)y (24) 

[SS]{t) = A^'(l-^+(^-2)yi(t)+l/2(t)), (25) 

[ISm) = N'(^-j^yi{t)+(^-^ + l^y2{t)-y3{t)y (26) 

The results above allows us to test if the equations for the moments are consistent with the 
pairwise ones. Starting from the the equation for the expected number of infecteds and upon 
using Eq. f|T8|) . we obtain. 



[i]it) = Nmit) = iV((/3 - 7)yi - y^) = N{(5{yi - y^) - 71/1) = ^[SI] - i[I]. (27) 

The calculations above can be repeated for all other equations and these confirm the one- 
to-one correspondence between the moment and pairwise equations. 
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3.3 Interpreting closures via moments 

The closure at the level of pairs and triples have to equally translate in a relation between 
the moments. First, the closure at the level of pairs is discussed. We look to use the pairwise 
closure [SI] = [S][I] = {N — [/])[/] to obtain a relation between the moments. The simple 
relation above in terms of the moments translates to 

N^y^it) - y,{t)) = N\y,it) - yl{t)) ^ y,{t) = yl{t). (28) 

For the triple closure the situation is different in that two different triples are closed, albeit 
using the same formal relation, and this could potentially lead to two different relations 
between the moments. The first triple closure, for the SSI triple, upon using Eqs. f l2T]) - fl23|) 
& fl25l) . leads to the following relation. 

The equation above, after some algebra, yields a closure at the level of the moments and 1/3 
can be given in function of the previous two moments as follows, 

1 ^ n ^ ^ iV - 2 (yi - y^f 

y. = -j^y. + (1 + j^)y. - • (so) 

It is worth noting that the closure for [ISI] yields the same closure relation for the third 
moment. 



4 The new, improved closure and its performance 

The novel closure put forward here is based on the empirical observation that Pkif) is well 
approximated by a binomial distribution B{n,p), where n and p depend on time and will 
specified in terms of the moments of the distribution. The first three moments of the binomial 
distribution can be specified easily in terms of the two parameters and are as follows, 

Yi = np (31) 
I2 = np + n{n — l)p^ (32) 
1^3 = np + ?)n[n — \)p'^ + n{n — l){n — 2)p^ . (33) 

Using Eqs. (^11) & (15^ . n and p can be expressed in term of Yi and Y2 as follows. 
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Plugging the expressions for p and n (Eq. [M]) into Eq. (1551) . the closure for the third moment 
is found to be 

Ys = -^-Y2- Y^{Y2 - Fi). (35) 

This relation defines the new triple closure and in terms of the density dependent moments 
is equivalent to 

?/3 = — + 4(^1 ~^2)- (36) 
yi N 

Using the equations for the first moment (|T8|) the closure at the level of the pairs yields the 
following approximate equation 

X, = - ^)xi - ^xl (37) 

Using the equations for the first two moments ( f[T5]) & flTI?]) ) and the closure at the level of 
the third moment yields 

xi = {(3 - j)xi - /3x2, (38) 
X2 = 2(/3-7)x2- 2/3x3 + ^((/3 + 7)a;i-/3a;2), (39) 

where ^ 

X3 = — - - xiX2 + —{xl - X2). (40) 
Xi JM 

Moreover, we can also define a simplified binomial closure by neglecting the order in 
the full binomial closure, provided also that xf — X2 is of 0{1/N). This will lead to 

X3 = 3xiX2 — 2x1, (41) 

with all the above in contrast with the the classic triple closure given by Eq. fl30p . 

The current setup allows us to compare the exact model as given by the Kolmogorov 
equations flKEp with transition rates given by Eq. ([1]) to three different approximate models. 
The first results from the pairwise closure and is given by Eq. ([5]). The second is a direct 
consequence of the closure at the level of triples and is given by Eqs. f|T0l) - f|T3l) . Finally, the 
third approximate system results from the novel binomial closure (see Eq. fl5^ ) with the 
approximate system defined by Eqs. fl55]) - fHn]) . The elegance of this approach stems from the 
fact that all numerical results are free from simulations and rely solely on the time integration 
of ODE systems. In Fig. 1, plots of the time evolution together with the approximation errors 
are given for the classic approximate models corresponding to pairwise and triple closure. 
The most significant feature of this plot is the order 1/A^ approximation error independently 
of the closure. The same approximation error suggests qualitative similarity between the 
two approximate models and it somewhat surprising that the closure at the level of triples 
produces no immediately obvious benefits or qualitatively different behaviour. In contrast. 
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in Fig. 2 the performance of the binomial closure is compared to the triple closure. The 
approximation error plot in Fig. 2b is the most significant as it shows numerical evidence 
that the binomial closure performs significantly better with error that is of order 1/iV^. 
While we were not yet able to prove analytically this result, for the case of pairwise closure 
we gave rigorous proof in [21 |TT], where the result is more general as apart from estimating 
the approximate error at the steady state also extends to the approximate error at different 
times. 

The quantitative comparison of the different closure relations is shown in Table 1. Her, 
the values of the steady states obtained for the different closures are compared to the exact 
value of the stationary prevalence. The advantage of comparing the steady state values is 
that these can be given analytically without solving ODEs numerically. The steady state 
from the pairwise closure Eq. is 



1 - 



7 
/3' 



(42) 



where (3 = tN. The steady state from the triple closure Eqs. (ITU]) - !!!^ is 



N ^ ' {N -1){N -2)f3 - N-f' 
The steady state from the binomial closure Eqs. fl38l) - fl40l) is 



(43) 



Nq' - 1 
N Nq-1'' 



(44) 



where q = 1 — 7//?. Finally, the value of the steady state from the Kolmogorov equations 
(KE) can be derived as 

N-l 



N-l 
k=0 



where Aq = 1 and 



7'=(A; + 1)A^*^ 
The approximation errors shown in Table 1 are 



N 



I ss 



N 



N 



N 
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In the Table these values are multiphed by 1000 to get_ numerical values closer to 1. Using the 
analytical expressions above it is easy to show that ^ ^ ^ and ^ ^ ^ when N ^ oo. 
However, we did not yet manage to formally prove that 





[I]ss 


/N, 


-I- ss 


[I]ss 


/N, 


I* 

ss 


[I]ss 


TV 


N 




N 


N 


N 



all tend to a non-zero number as N gets large and thus analytically confirming the 0{1/N) 
and 0{1/N'^) errors. 



5 Discussion 

In this paper, we proposed a novel/new closure that produces qualitatively different results 
when compared to the classic closure. Namely, the difference between the exact system 
from the solution of the approximate model that results from the binomial closure is smaller 
(i.e. 0{1/N'^) compared to 0{1/N)) and this has been illustrated via numerical examples. 
The paper has also identified the link between the equations of the moments and pairwise 
equations. In particular, it is worth noting that in terms of moments, the pairwise equations 
reduce to two simple ODEs and this is in contrast with four equations for the pairwise model. 
However, this apparent discrepancy is not as marked since the pairwise model can be easily 
reduced to three equations by noting that, for example, [SS] = {N — 1){N — [I]) — [SI]. This 
still leaves one extra equation in the pairwise model when compared to the equivalent system 
in terms of moments. This suggests that the three equations are not completely independent 
in some sense that is not easy to define or pinpoint. 

The present study also allows us to make a few more important observations. For ex- 
ample, one can attempt to close the moment equations at the second moment but, without 
using the equivalent closure that follows from the pairwise model. This can be achieved by 
considering the limiting case of n = N and expressing the second moment in terms of the 
first. This will lead to a new closure at the level of pairs but, this performs similarly to the 
classic closure at the pair level. 

The procedure that we illustrate in this paper can be generalised and applied within 
the general context of model reduction, where Kolmogorov-typc systems arc reduced to 
approximate models with far fewer equations. Our method is based on a semi-heuristic ar- 
gument which relies on the assumption that the distribution of the stochastic process in time 
can be well-approximated by some theoretical distribution with time-dependent parameters. 
Moreover, the distribution has to be such that there is a unique, one-to-one correspondence 
between the moments and the parameters of the distribution. If the conditions above hold, 
the theoretical distribution will determine the exact type of closure that is necessary to ob- 
tain a reduction in the number of equations. The method also allow us to close at different 
levels or moments and to explore the benefits of closing at higher order moments. While 
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the initial results and applicability of our method is promising (see example below for the 
homogeneous random graph case), we are yet to fully explore the applicability of our method 
or to identify the sufficient and/or necessary requirements that guarantees that the proposed 
closure or reduction procedure will be applicable. It is worth noting that here the main focus 
in not on the approximation of simulation results corresponding to a stochastic process but 
rather the model reduction where the Kolmogorov equations are already given. 

Our result about the new closure relation applies directly to homogeneous random graphs 
provided that the Kolmogorov equations for this case take the form given in Eq. (KE) with 
coefficients defined in Eq. ([2]). The formulation of the result is even more straightforward if 
the following slight modification for the transmission rate is introduced 

N -k 
ttk = rnk ^ . 

In this case, the coefficients and in Eq. ([T]) and Eq. are formally equivalent, i.e. 
both can be written as = -^k{N — k) with 5 = l3 for the complete graph, and 6 = rn for 
the homogeneous random graph. It is worth noting that, in the case of the complete graph, 
r is of order 1/N and /3 is independent of A^, while for the homogeneous random graph, r 
is independent of A^. However, in both cases the coefficient can be given as the product 
of k{N — k) with a term. Therefore in the case of a homogeneous random graph the 
equation for the moments can be easily obtained from Eqs. flTSl) & f[T^ by writing rn instead 
of p. Thus the ODE system with the binomial closure for the homogeneous random graph 
is the same as Eqs. (l38|) - fH0|) following a change that again replaces (3 with rn. 

While there are some limitations to the straightforward extension of our results to arbi- 
trary graphs and/or dynamics, the present paper provides the basis for some more rigorous 
analysis as well as the possibility of deriving new closures. As such, this study offers a new or 
different direction that can lead to further progress in the area of approximating stochastic 
processes independently of whether these involve graphs. 
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N 


100 


200 


400 


800 


1000 X pair appr. err. 


6.9486 


3.4008 


1.6832 


0.8374 


1000 X triple appr. err. 


1.2355 


0.5729 


0.2763 


0.1357 


1000 X binomial appr. err. 


0.1689 


0.0395 


0.0096 


0.0024 



Table 1: Numerical values for the discrepancy between the exact and approximate solutions 
for different network sizes. 
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Figure 1: (a) Time evolution of the fraction infected (prevalence) based on the exact model 
and two different closures; pair and triple level. Here, N = 200, 7 = 2 and /3 = 5. (b) 
Approximation error (absolute value of the difference between the exact and two approximate 
models at steady state) plotted for different system sizes for the same transmission and 
recovery parameter values as at (a). 
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Figure 2: (a) Time evolution of the fraction infected (prevalence) based on the exact model 
and three different closures; pair, triple level and the binomial closure. Here, iV = 200, 7 = 2 
and /3 = 5. (b) Approximation error (absolute value of the difference between the exact and 
approximate models at steady state given by the triple and binomial closures) plotted for 
different system sizes for the same transmission and recovery parameter values as at (a). It 
is worth noting that the x-axis in (b) is in terms of as opposed to as in Fig. 1(b). 
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